rm(list=ls())

in.dir = "/home1/99/jc152199/In/"
out.dir = "/home1/99/jc152199/Out/"

#Set working directory

setwd(in.dir)

#load necessary libraries

library("SDMTools")
library("quantreg")
library("adehabitat")

#Import the dataset for analysis

samples= read.csv("samplesmaster.csv",header=T)
georef = read.csv("georefmaster.csv", header=T)
s.fields = names(samples)
g.fields = names(georef)

#Perform quantile regression for each of 4 species against it's own ES values

cr.rq=rq(georef$cr_avgabund~georef$cr_es,data=georef,tau=.95,method="br",model=T)
sb.rq=rq(georef$sb_avgabund~georef$sb_es+0,data=georef,tau=.975,method="br",model=T)
gq.rq=rq(georef$gq_avgabund~georef$gq_es+0,data=georef,tau=.975,method="br",model=T)
lc.rq=rq(georef$lc_avgabund~georef$lc_es,data=georef,tau=.975,method="br",model=T)

#Make a data frame from the regression residuals and the raw georef data

master=data.frame(georef, cr.rq$residuals, gq.rq$residuals,lc.rq$residuals, sb.rq$residuals)
#master$sbmax = ((master$sb_es*sb.rq$coefficients[1]))
#master2= subset(master, master$lsc_mean>0 & master$LAI4_forest>0)


#Create the multiple linear model to explain residual values

gq.lm = lm(master$gq.rq.residuals~master$MODIS_mean, data=master)
lc3.lm = lm(master$lc.rq.residuals~master$bc05micro, data=master)
sb1.lm = lm(master2$sb.rq.residuals~master2$lsc_mean+master$LAI4_forest, data=master)
sb2.lm = lm(master2$sb.rq.residuals~master2$lsc_mean+master$LAI4_forest+master$lc_es, data=master)

#Read in spatial data

sb.asc = read.asc(paste(in.dir,"SAPBASI.asc",sep=""))
lc.asc = read.asc(paste(in.dir,"LAMCOGG.asc",sep=""))
bc.asc = read.asc(paste(in.dir,"bc_05_buff.asc",sep=""))

#Perform functions on two linear models to get spatial residual and max abundance values

#master2$sbresid = (sb1.lm$coefficients[1]+(sb1.lm$coefficients[2]*master2$lsc_mean)+(sb1.lm$coefficients[3]*master2$LAI4_forest))
#master3$sbresid = (sb2.lm$coefficients[1]+(sb2.lm$coefficients[2]*master3$lsc_mean)+(sb2.lm$coefficients[3]*master3$LAI4_forest)+(sb2.lm$coefficients[4]*master3$lc_es))
#master$sbmax = ((master$sb_es*sb.rq$coefficients[1]))
#master2$sbreal = master2$sbmax-master2$sbresid

master$gqresid = (gq.lm$coefficients[1]+(gq.lm$coefficients[2]*master$MODIS_mean)
master$gqmax = (gq.rq$coefficients[1]*master$gq_es)
master$gqreal = (master$gqmax-master$gqresid)

summary(master$gqresid)
summary(master$gqmax)
summary(master$gqreal)

#Column bind max and real data to master to create new scatterplots and perform new OLS and quantregs
#First subset master to get just XY data

#XY = data.frame(master$east,master$north)

#max = extract.data(XY, sbmax)
#real = extract.data(XY, sbreal)

#Now perform scatterplots of master$lc_es against max and real

png(paste(out.dir,"realabund1.png", sep=""), height=800, width=800)

plot(master$sb_es,master$sbreal1, xlab="ES", ylab="RealAbund1", main="SAPBASI")

dev.off()

png(paste(out.dir,"maxabund.png", sep=""), height=800, width=800)

plot(master$sb_es,sbmax, xlab="ES", ylab="MaxAbund", main="SAPBASI")

dev.off()

png(paste(out.dir,"realabund2.png", sep=""), height=800, width=800)

plot(master$sb_es,sbreal2, xlab="ES", ylab="realabund2", main="SAPBASI")

dev.off()




